global main "C:/Users/silvio/Documents/CVR/reply_to_mvb"
global excel "$main/excel"  
global code "$main/code"  
global data "$main/data"  

global datafile "datacvr" 
global name "ballman_mijk_allstrata"

matrix drop _all
set more off
set matsize 1392
set seed 98034

cd $code

//User defined macros end

use "${data}/${datafile}" , replace
mkmat y


**************************************************************
matrix am=J(8,1,0)
matrix one7=J(1,7,1)
matrix one8=J(1,8,1)
matrix one58=J(1,58,1)
matrix truenobs=J(1,3,0)
**************************************************************

matrix m_matsipr_all=J(8,7,0)

matrix freqneg=J(9,2,0)
*********************************
//Strata
*********************************
forval ii = 1/58{  /* begin loop strata */
scalar jhat=`ii'
di jhat

matrix m_mate=J(8,1,0)
matrix m_mats=J(8,1,0)
matrix m_matsi=J(8,1,0)

matrix m_mate0=J(8,1,0)
matrix m_mats0=J(8,1,0)

matrix m_matsi0=J(8,1,0)

scalar d_mats=0
scalar d_matsi=0

scalar p_mats=0
scalar p_matsi=0

matrix m_matsi_all=J(8,7,0)




*******************************
//nobs, N direct, indirect
*******************************
scalar perpeh=1
scalar jj1=(jhat-1)*8*3+(perpeh-1)*8+1
scalar jj2=jj1+6
matrix a=y[jj1..jj2,1]              /*original data*/
run "$code/analya"
matrix mijk0= mijk
***************************************************************
//Integers for mijk
forval ij = 1/7 {       //models
forval ik = 1/8 {      //Counts
matrix mijk[`ik',`ij']=round(mijk[`ik',`ij']) 
}
}
***************************************************************

scalar modeli_e=modeli
scalar imodeli_e=1
scalar chi2_e=0
if modeli>0{
scalar imodeli_e=0
matrix m_mate=mijk[1..8,modeli]           /*estimated parameters*/
matrix m_mate0=mijk0[1..8,modeli]           /*estimated parameters*/
}
matrix ae=a
scalar modeli_e=modeli






scalar perpeh=2
scalar jj1=(jhat-1)*8*3+(perpeh-1)*8+1
scalar jj2=jj1+6
matrix a=y[jj1..jj2,1]              /*original data*/
run "$code/analya"
matrix mijk0= mijk

***************************************************************
//Integers for mijk
forval ij = 1/7 {       //models
forval ik = 1/8 {      //Counts
matrix mijk[`ik',`ij']=round(mijk[`ik',`ij']) 
}
}
***************************************************************

scalar chi2_s=0
scalar modeli_s=modeli
scalar imodeli_s=1
if modeli>0{
scalar imodeli_s=0
matrix m_mats=mijk[1..8,modeli]           /*estimated parameters*/
matrix m_mats0=mijk0[1..8,modeli]           /*estimated parameters*/

scalar chi2_s=chi2ml[modeli,1]
scalar d_mats=dfml[modeli,1]
scalar p_mats=prm[modeli,1]

}
matrix as=a


scalar iperf_mats0=0
scalar iperf_mats=0
forval i = 1/7{
if a[`i',1]==m_mats0[`i',1]{
scalar iperf_mats0=iperf_mats0+1
}
if a[`i',1]==m_mats[`i',1]{
scalar iperf_mats=iperf_mats+1
}
}
if p_mats ~=. & p_mats >=prlow & p_mats <= prhigh {
scalar ip_mats=0
}
else {
scalar ip_mats=1
}



*E+S e
matrix a=ae+as                   /*5. Simulated data of sum of E+S */
matrix mat_norm=m_mats
run "$code/analya"
matrix mijk0= mijk


***************************************************************
//Integers for mijk
forval ij = 1/7 {       //models
forval ik = 1/8 {      //Counts
matrix mijk[`ik',`ij']=round(mijk[`ik',`ij']) 
}
}
***************************************************************

scalar modeli_ei=modeli
scalar imodeli_ie=1
if modeli>0{
scalar imodeli_ie=0
matrix m_mates=mijk[1..8,modeli]           /*estimated parameters*/
matrix m_mates0=mijk0[1..8,modeli]           /*estimated parameters*/

}
*E+S s
matrix a=ae+as                   /*5. Simulated data of sum of E+S */
matrix mat_norm=m_mate
run "$code/analya"
matrix mijk0= mijk


***************************************************************
//Integers for mijk
forval ij = 1/7 {       //models
forval ik = 1/8 {      //Counts
matrix mijk[`ik',`ij']=round(mijk[`ik',`ij']) 
}
}
***************************************************************

scalar modeli_si=modeli
scalar imodeli_is=1
if modeli>0{
scalar imodeli_is=0
matrix m_mates=mijk[1..8,modeli]           /*estimated parameters*/
if modeli_e>0{
scalar d_matsi=dfml[modeli,1]
matrix m_matsi=m_mates-m_mate                /*indirect parameters*/
matrix m_matsi0=m_mates0-m_mate0                /*indirect parameters*/
matrix m_matsi_all=mijk-m_mate*one7
}
}
matrix m_mate7=m_mate[1..7,1] 
matrix m_mats7=m_mats[1..7,1] 

matrix m_matsi7=m_matsi[1..7,1] 
*******************************
//nobs, N direct, indirect end
*******************************
//Direct observations for E,S and ES
//Observations
matrix truenobs[1,1]=one7*as
//Totals Direct
matrix truenobs[1,2]=one8*m_mats 
//Totals Indirect
matrix truenobs[1,3]=one8*m_matsi 

scalar indinegs=0
forval i = 1/8{
if m_matsi[`i',1]<0{
scalar indinegs=1
}
}

scalar indinegs_at=0
forval j = 1/7{
scalar ixxs=0
forval i = 1/8{
if m_matsi_all[`i',`j']<0{
scalar ixxs=1
matrix m_matsipr_all[`i',`j']=m_matsipr_all[`i',`j']+1
}
}
scalar indinegs_at=indinegs_at+ixxs
}

scalar indinegs_ats=int(indinegs_at/7)


scalar iperf_matsi0=0
matrix mijk=m_matsi0[1..7,1]
matrix a=as
scalar chi2=0
forval i = 1/7{

if a[`i',1]==mijk[`i',1]{
scalar iperf_matsi0=iperf_matsi0+1
}


if a[`i',1]~=0 | mijk[`i',1]~=0{
scalar chi2=chi2+(a[`i',1]-mijk[`i',1])*(a[`i',1]-mijk[`i',1])/mijk[`i',1]
}
} /* end loop */

scalar chi2_is=chi2
scalar p_matsi=1 - chi2(d_matsi,chi2_is)
scalar c_mats=chi2_s/d_mats
scalar c_matsi=chi2_is/d_matsi

if p_matsi ~=. & p_matsi >=prlow & p_matsi <= prhigh {
scalar ip_matsi=0
}
else {
scalar ip_matsi=1
}


scalar iperf_mats0=int(iperf_mats0/7)
scalar iperf_matsi0=int(iperf_matsi0/7)
scalar indinegs_perf0=max(indinegs,iperf_matsi0)
scalar indinegs_perf0_p=max(ip_matsi,indinegs_perf0)

if jhat==1{
matrix res_s=jhat,2,truenobs[1,1..3],as',m_mats0',modeli_s,m_matsi0',modeli_si,chi2_s,d_mats,c_mats,p_mats,chi2_is,d_matsi,c_matsi,p_matsi,indinegs_ats, indinegs, iperf_matsi0, indinegs_perf0, ip_matsi, indinegs_perf0_p
}
else{
matrix res_s=res_s\jhat,2,truenobs[1,1..3],as',m_mats0',modeli_s,m_matsi0',modeli_si,chi2_s,d_mats,c_mats,p_mats,chi2_is,d_matsi,c_matsi,p_matsi,indinegs_ats, indinegs, iperf_matsi0, indinegs_perf0, ip_matsi, indinegs_perf0_p
}
} //Strata


clear
svmat res_s, names(v)
ren v1 stratum		
ren v2 perpe
ren v3 nobs
ren v4 nobstot
ren v5 nobstoti

ren v6 m100
ren v7 m010
ren v8 m110
ren v9 m001
ren v10 m101
ren v11 m011
ren v12 m111

ren v13 d100
ren v14 d010
ren v15 d110
ren v16 d001
ren v17 d101
ren v18 d011
ren v19 d111
ren v20 d000

ren v21 modeli

ren v22 i100
ren v23 i010
ren v24 i110
ren v25 i001
ren v26 i101
ren v27 i011
ren v28 i111
ren v29 i000

ren v30 modeli_i

ren v31 chi_d
ren v32 d_d
ren v33 chidf_d
ren v34 p_d

ren v35 chi_i
ren v36 d_i
ren v37 chidf_i
ren v38 p_i

ren v39 indinegs_ats
ren v40 indineg
ren v41 iper0_i
ren v42 indineg_perf0
ren v43 ip_i
ren v44 indineg_perf0_p

export excel * ///
using "$excel/${name}.xls", sheet("Table 3 and 4. Main Text") sheetmodify firstrow(varlabels) missing(".")
clear

matrix m_matsi_sint=res_s[1..58,39..44]
matrix m_matsi_sint=one58*m_matsi_sint

svmat m_matsi_sint, names(v)
ren v1 indinegs_ats
ren v2 indineg
ren v3 iper0_i
ren v4 indineg_perf0
ren v5 ip_i
ren v6 indineg_perf0_p

export excel * ///
using "$excel/${name}.xls", sheet("Table 5 Main Text") sheetmodify firstrow(varlabels) missing(".")
clear

svmat m_matsipr_all, names(s)
export excel * ///
using "$excel/${name}.xls", sheet("Table 5 Supplemental Material") sheetmodify firstrow(varlabels) missing(".")
clear


